Model Report
Sam Abbott
2018-09-07
Load packages required for analysis
Load the fitted model and prior draws.
model <- ModelTBBCGEngland::read_libbi(file.path(model_dir, "libbi", "posterior"))
priors <- readRDS(file.path(model_dir, "data", "prior-params.rds"))
model## Wrapper around LibBi
## ======================
## Model: Baseline
## Run time: 3790.562 seconds
## Number of samples: 1000
## State trajectories recorded: E foi H L N P S T_E T_P YearlyAgeCases YearlyCases YearlyEAgeCases YearlyECases YearlyEPulCases YearlyEPulDeaths YearlyNonUKborn YearlyPAgeCases YearlyPCases YearlyPulCases YearlyPulDeaths
## Observation trajectories recorded: YearlyAgeInc YearlyHistPInc YearlyInc
## Parameters recorded: alpha_t c_eff c_hist chi_init delta epsilon_h_0_4 epsilon_h_15_89 epsilon_h_5_14 epsilon_l_0_4 epsilon_l_15_89 epsilon_l_5_14 HistMeasError kappa_0_4 kappa_15_89 kappa_5_14 M mu_e_0_14 mu_e_15_59 mu_e_60_89 mu_p_0_14 mu_p_15_59 mu_p_60_89 nu_e_0_14 nu_e_15_89 nu_p_0_14 nu_p_15_89 phi_0_14 phi_15_59 phi_60_89 rho_0_14 rho_15_59 rho_60_89 Upsilon_0_14 Upsilon_15_59 Upsilon_60_89 zeta_e_0_14 zeta_e_15_59 zeta_e_60_89 zeta_p_0_14 zeta_p_15_59 zeta_p_60_89- Get the rejection rates of parameters
traces <- coda::mcmc(rbi::get_traces(model))
coda::rejectionRate(traces)## M c_eff c_hist delta
## 0.989990 0.989990 0.989990 0.989990
## epsilon_h_0_4 epsilon_h_5_14 epsilon_h_15_89 kappa_0_4
## 0.989990 0.989990 0.989990 0.989990
## kappa_5_14 kappa_15_89 epsilon_l_0_4 epsilon_l_5_14
## 0.989990 0.989990 0.989990 0.989990
## epsilon_l_15_89 phi_0_14 phi_15_59 phi_60_89
## 0.989990 0.989990 0.989990 0.989990
## Upsilon_0_14 Upsilon_15_59 Upsilon_60_89 rho_0_14
## 0.990991 0.989990 0.990991 0.989990
## rho_15_59 rho_60_89 nu_p_0_14 nu_p_15_89
## 0.989990 0.989990 0.989990 0.989990
## nu_e_0_14 nu_e_15_89 zeta_p_0_14 zeta_p_15_59
## 0.989990 0.989990 0.989990 0.989990
## zeta_p_60_89 zeta_e_0_14 zeta_e_15_59 zeta_e_60_89
## 0.989990 0.989990 0.989990 0.989990
## mu_p_0_14 mu_p_15_59 mu_p_60_89 mu_e_0_14
## 0.989990 0.989990 0.989990 0.989990
## mu_e_15_59 mu_e_60_89 chi_init alpha_t.0
## 0.989990 0.989990 0.989990 1.000000
## alpha_t.1 alpha_t.2 alpha_t.3 alpha_t.4
## 1.000000 1.000000 1.000000 1.000000
## alpha_t.5 HistMeasError
## 1.000000 0.989990- Evaluate Posterior Traces
plot_param(model, scales = "free", plot_type = "trace") + theme(legend.position = "none")## Aggregate function missing, defaulting to 'length'
- Plot overview of prior and posterior densities
plot_param(model, prior_params = priors, scales = "free")## Aggregate function missing, defaulting to 'length'
- Table of priors and posteriors
param_sum <- plot_param(model, prior_params = priors, plot_data = FALSE) %>%
group_by(Distribution, parameter, length) %>%
summarise(median = median(value),
lll = quantile(value, 0.025),
hhh = quantile(value, 0.975)) %>%
mutate(value = pretty_ci(median, lll, hhh, sep = ", ")) %>%
group_by(Distribution, parameter) %>%
mutate(Parameter = case_when(max(length) > 1 ~ paste(parameter, 1:n(), sep = "-"),
TRUE ~ as.character(parameter))
) %>%
ungroup %>%
select(Distribution, Parameter, value) %>%
spread(key = "Distribution", value = "value") %>%
select(Parameter, Prior, Posterior)## Aggregate function missing, defaulting to 'length'saveRDS(param_sum, file.path(model_dir, "data", "sum_prior_posterior.rds"))
knitr::kable(param_sum)| Parameter | Prior | Posterior |
|---|---|---|
| alpha_t-1 | 0.85 (0.77, 0.90) | 0.86 (0.86, 0.86) |
| alpha_t-2 | 0.69 (0.52, 0.80) | 0.75 (0.75, 0.75) |
| alpha_t-3 | 0.57 (0.33, 0.71) | 0.45 (0.45, 0.45) |
| alpha_t-4 | 0.57 (0.37, 0.71) | 0.54 (0.54, 0.54) |
| alpha_t-5 | 0.25 (-0.11, 0.48) | 0.09 (0.09, 0.09) |
| alpha_t-6 | 0.21 (-0.44, 0.55) | 0.14 (0.14, 0.14) |
| c_eff | 2.61 (0.14, 4.85) | 4.87 (4.87, 4.87) |
| c_hist | 12.48 (10.16, 14.87) | 10.82 (10.82, 10.82) |
| chi_init | 0.19 (0.08, 0.28) | 0.25 (0.25, 0.25) |
| delta | 0.78 (0.70, 0.87) | 0.70 (0.70, 0.70) |
| epsilon_h_0_4 | 1.49 (0.09, 5.16) | 4.40 (4.40, 4.40) |
| epsilon_h_15_89 | 24.65 (1.42, 73.48) | 94.34 (94.34, 94.38) |
| epsilon_h_5_14 | 3.90 (0.18, 11.83) | 9.50 (9.50, 9.51) |
| epsilon_l_0_4 | 616.42 (33.73, 1720.16) | 51.48 (51.41, 51.52) |
| epsilon_l_15_89 | 1081.98 (58.39, 3378.35) | 1766.49 (1765.83, 1766.49) |
| epsilon_l_5_14 | 489.66 (28.10, 1488.47) | 921.11 (921.11, 922.35) |
| HistMeasError | 1.00 (0.61, 1.40) | 1.25 (1.25, 1.25) |
| kappa_0_4 | 0.86 (0.05, 2.79) | 0.67 (0.67, 0.67) |
| kappa_15_89 | 1.09 (0.05, 3.41) | 2.90 (2.90, 2.90) |
| kappa_5_14 | 0.97 (0.06, 3.10) | 1.44 (1.44, 1.44) |
| M | 0.24 (0.01, 0.49) | 0.49 (0.49, 0.49) |
| mu_e_0_14 | 275.68 (217.01, 341.37) | 239.52 (239.52, 239.52) |
| mu_e_15_59 | 193.36 (72.76, 318.29) | 401.85 (401.72, 401.85) |
| mu_e_60_89 | 37.15 (2.08, 103.08) | 71.24 (71.24, 71.31) |
| mu_p_0_14 | 241.67 (155.25, 324.07) | 190.08 (190.07, 190.08) |
| mu_p_15_59 | 82.93 (4.33, 264.25) | 248.62 (248.56, 248.62) |
| mu_p_60_89 | 50.54 (3.38, 158.42) | 39.28 (39.27, 39.28) |
| nu_e_0_14 | 0.19 (0.00, 1.16) | 0.19 (0.19, 0.19) |
| nu_e_15_89 | 0.31 (0.01, 1.93) | 1.15 (1.15, 1.15) |
| nu_p_0_14 | 0.11 (0.00, 0.78) | 0.01 (0.01, 0.01) |
| nu_p_15_89 | 0.24 (0.01, 1.25) | 4.40 (4.40, 4.40) |
| phi_0_14 | 0.57 (0.28, 1.04) | 0.53 (0.53, 0.53) |
| phi_15_59 | 0.61 (0.26, 1.18) | 1.72 (1.72, 1.72) |
| phi_60_89 | 0.59 (0.26, 1.11) | 0.50 (0.50, 0.50) |
| rho_0_14 | 0.30 (0.26, 0.34) | 0.27 (0.27, 0.27) |
| rho_15_59 | 0.65 (0.64, 0.66) | 0.69 (0.69, 0.69) |
| rho_60_89 | 0.54 (0.52, 0.55) | 0.56 (0.56, 0.56) |
| Upsilon_0_14 | 0.63 (0.61, 0.65) | 0.60 (0.60, 0.60) |
| Upsilon_15_59 | 0.71 (0.70, 0.71) | 0.70 (0.70, 0.70) |
| Upsilon_60_89 | 0.75 (0.74, 0.76) | 0.76 (0.76, 0.76) |
| zeta_e_0_14 | 122.39 (57.91, 187.56) | 142.61 (142.56, 142.61) |
| zeta_e_15_59 | 60.81 (3.67, 171.01) | 41.51 (41.51, 41.58) |
| zeta_e_60_89 | 135.73 (55.90, 212.67) | 98.79 (98.78, 98.85) |
| zeta_p_0_14 | 94.95 (18.70, 179.68) | 19.81 (19.81, 19.89) |
| zeta_p_15_59 | 78.52 (3.22, 242.86) | 605.28 (605.23, 605.35) |
| zeta_p_60_89 | 121.96 (15.39, 248.56) | 63.16 (63.13, 63.16) |
- Plot an overview of states
plot_state(model, summarise = TRUE)
- Plot an overview of states zoomed into 1990 onwards
plot_state(model, summarise = TRUE, start_time = 59)
- Plot fitted observed states.
plot_state(model, states = c("YearlyHistPInc", "YearlyInc"), summarise = TRUE, start_time = 59)
- Table fitted states
pred_obs <- plot_state(model, states = c("YearlyHistPInc", "YearlyInc", "YearlyAgeInc"), summarise = FALSE, start_time = 59, plot_data = FALSE)
obs <- ModelTBBCGEngland::setup_model_obs() %>%
bind_rows(.id = "state")
pred_obs <- pred_obs %>%
dplyr::filter(Average == "median") %>%
mutate(pred_value = pretty_ci(Count, lll, hhh, digits = 0)) %>%
left_join(obs, by = c("state" = "state", "time" = "time", "age" = "age")) %>%
select(Observation = state, Age = age, Year = time, `Observed Incidence` = value, `Predicted Incidence` = pred_value) %>%
mutate(Year = Year + 1931)
sum_obs_table <- pred_obs %>%
dplyr::filter(Observation %in% c("YearlyInc", "YearlyHistPInc")) %>%
drop_na(`Observed Incidence`) %>%
select(-Age) %>%
mutate(Observation = case_when(Observation == "YearlyInc" ~ "All TB cases",
Observation == "YearlyHistPInc" ~ "Pulmonary TB cases"))
saveRDS(sum_obs_table, file.path(model_dir, "data", "sum_obs_table.rds"))
age_cases_table <- pred_obs %>%
dplyr::filter(Observation %in% c("YearlyAgeInc")) %>%
drop_na(`Observed Incidence`) %>%
select(-Observation) %>%
group_by(Year) %>%
mutate(Age = c(paste0(seq(0, 65, 5), "-", seq(4, 69, 5)), "70-89")) %>%
ungroup
saveRDS(age_cases_table, file.path(model_dir, "data", "age_cases_table.rds"))knitr::kable(sum_obs_table)| Observation | Year | Observed Incidence | Predicted Incidence |
|---|---|---|---|
| Pulmonary TB cases | 1990 | 3618 | 4154 (4147 to 4366) |
| Pulmonary TB cases | 1991 | 3596 | 4063 (4062 to 4106) |
| Pulmonary TB cases | 1992 | 3816 | 3804 (3804 to 4009) |
| Pulmonary TB cases | 1993 | 3961 | 3748 (3748 to 3877) |
| Pulmonary TB cases | 1994 | 3694 | 3657 (3608 to 3806) |
| Pulmonary TB cases | 1995 | 3711 | 3608 (3595 to 3723) |
| Pulmonary TB cases | 1996 | 3690 | 3581 (3544 to 3700) |
| Pulmonary TB cases | 1997 | 3947 | 3632 (3566 to 3717) |
| Pulmonary TB cases | 1998 | 3926 | 3466 (3466 to 3594) |
| Pulmonary TB cases | 1999 | 3827 | 3506 (3476 to 3578) |
| All TB cases | 2000 | 1803 | 1811 (1774 to 1937) |
| All TB cases | 2001 | 1866 | 1817 (1766 to 1817) |
| All TB cases | 2002 | 1833 | 1779 (1684 to 1784) |
| All TB cases | 2003 | 1685 | 1652 (1612 to 1772) |
| All TB cases | 2004 | 1776 | 1688 (1669 to 1749) |
knitr::kable(age_cases_table)| Age | Year | Observed Incidence | Predicted Incidence |
|---|---|---|---|
| 0-4 | 2000 | 75 | 95 (79 to 98) |
| 5-9 | 2000 | 59 | 55 (55 to 74) |
| 10-14 | 2000 | 75 | 74 (64 to 77) |
| 15-19 | 2000 | 112 | 14 (14 to 23) |
| 20-24 | 2000 | 133 | 35 (26 to 36) |
| 25-29 | 2000 | 116 | 60 (54 to 66) |
| 30-34 | 2000 | 147 | 75 (60 to 90) |
| 35-39 | 2000 | 99 | 121 (102 to 121) |
| 40-44 | 2000 | 83 | 146 (143 to 166) |
| 45-49 | 2000 | 105 | 192 (155 to 192) |
| 50-54 | 2000 | 116 | 204 (189 to 204) |
| 55-59 | 2000 | 112 | 228 (220 to 228) |
| 60-64 | 2000 | 101 | 228 (209 to 236) |
| 65-69 | 2000 | 99 | 182 (182 to 214) |
| 70-89 | 2000 | 371 | 137 (137 to 145) |
| 0-4 | 2001 | 113 | 84 (78 to 84) |
| 5-9 | 2001 | 65 | 52 (52 to 67) |
| 10-14 | 2001 | 51 | 82 (70 to 89) |
| 15-19 | 2001 | 130 | 19 (13 to 24) |
| 20-24 | 2001 | 145 | 41 (34 to 44) |
| 25-29 | 2001 | 127 | 46 (46 to 72) |
| 30-34 | 2001 | 122 | 71 (71 to 92) |
| 35-39 | 2001 | 128 | 111 (96 to 127) |
| 40-44 | 2001 | 107 | 152 (128 to 164) |
| 45-49 | 2001 | 101 | 185 (175 to 194) |
| 50-54 | 2001 | 96 | 213 (198 to 221) |
| 55-59 | 2001 | 91 | 232 (206 to 233) |
| 60-64 | 2001 | 94 | 214 (184 to 224) |
| 65-69 | 2001 | 105 | 198 (196 to 235) |
| 70-89 | 2001 | 391 | 151 (121 to 151) |
| 0-4 | 2002 | 102 | 85 (69 to 100) |
| 5-9 | 2002 | 62 | 67 (54 to 71) |
| 10-14 | 2002 | 64 | 75 (66 to 75) |
| 15-19 | 2002 | 129 | 16 (12 to 17) |
| 20-24 | 2002 | 148 | 37 (25 to 37) |
| 25-29 | 2002 | 120 | 46 (44 to 55) |
| 30-34 | 2002 | 112 | 76 (64 to 83) |
| 35-39 | 2002 | 127 | 86 (86 to 115) |
| 40-44 | 2002 | 98 | 138 (127 to 151) |
| 45-49 | 2002 | 89 | 175 (159 to 176) |
| 50-54 | 2002 | 104 | 189 (189 to 219) |
| 55-59 | 2002 | 116 | 203 (203 to 217) |
| 60-64 | 2002 | 88 | 215 (200 to 223) |
| 65-69 | 2002 | 90 | 195 (190 to 200) |
| 70-89 | 2002 | 384 | 145 (131 to 159) |
| 0-4 | 2003 | 74 | 102 (80 to 102) |
| 5-9 | 2003 | 46 | 51 (51 to 71) |
| 10-14 | 2003 | 59 | 55 (55 to 72) |
| 15-19 | 2003 | 102 | 15 (11 to 17) |
| 20-24 | 2003 | 116 | 24 (24 to 28) |
| 25-29 | 2003 | 137 | 54 (49 to 54) |
| 30-34 | 2003 | 118 | 66 (66 to 73) |
| 35-39 | 2003 | 113 | 108 (103 to 108) |
| 40-44 | 2003 | 109 | 115 (115 to 149) |
| 45-49 | 2003 | 83 | 185 (153 to 185) |
| 50-54 | 2003 | 102 | 188 (184 to 200) |
| 55-59 | 2003 | 92 | 203 (190 to 212) |
| 60-64 | 2003 | 98 | 203 (184 to 208) |
| 65-69 | 2003 | 94 | 224 (182 to 224) |
| 70-89 | 2003 | 342 | 126 (126 to 145) |
| 0-4 | 2004 | 112 | 106 (68 to 106) |
| 5-9 | 2004 | 71 | 73 (57 to 76) |
| 10-14 | 2004 | 81 | 98 (72 to 98) |
| 15-19 | 2004 | 111 | 19 (14 to 19) |
| 20-24 | 2004 | 120 | 35 (11 to 35) |
| 25-29 | 2004 | 136 | 37 (37 to 48) |
| 30-34 | 2004 | 129 | 45 (45 to 71) |
| 35-39 | 2004 | 103 | 99 (87 to 99) |
| 40-44 | 2004 | 128 | 120 (110 to 137) |
| 45-49 | 2004 | 94 | 157 (148 to 186) |
| 50-54 | 2004 | 103 | 177 (177 to 196) |
| 55-59 | 2004 | 98 | 221 (195 to 221) |
| 60-64 | 2004 | 86 | 214 (193 to 214) |
| 65-69 | 2004 | 92 | 194 (192 to 202) |
| 70-89 | 2004 | 312 | 150 (131 to 150) |
- Plot age distributed cases.
plot_state(model, states = c("YearlyAgeInc"), start_time = 59) + theme(legend.position = "none")